Physiology Analysis

Script to analyze larval size, symbiont density, and examine correlations between physiological responses.

Setup

Set up workspace, set options, and load required packages.

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

1. Larval Size

Import and manipulate data

# Larval size data
size <- read_csv("Mcap2020/Data/Physiology/Size/larval_size.csv")

#metadata
metadata <- read_csv("Mcap2020/Data/lifestage_metadata.csv")

size <- left_join(size, metadata)
size$hpf <- as.factor(size$hpf)

Prep data frame.

# Calculate mean counts for each sample
size <- size %>%
  dplyr::select(tube.ID, lifestage, replicate, `area (mm)`, hpf, group)%>%
  drop_na()%>% #remove na's that could not be measured
  rename(area=`area (mm)`) #rename column

size$lifestage<-as.factor(size$lifestage)

Plotting

Plot data with mean and standard error for each lifestage.

size %>%
  ggplot(aes(x = lifestage, y = area, color = lifestage)) +
  labs(x = "",y = "Mean Larval Size (mm^2)") +
geom_jitter(width = 0.1) +                                            # Plot all points
  stat_summary(fun.data = mean_cl_normal, fun.args = list(mult = 1),    # Plot standard error
               geom = "errorbar", color = "black", width = 0.5) +
  stat_summary(fun = mean, geom = "point", color = "black") + # Plot mean
  theme_classic()

Present means and standard error of each group and save summary table

size%>%
  group_by(lifestage, hpf)%>%
  summarise(n=length(area),
            Mean=format(round(mean(area), 3), 3), 
            SE=format(round(sd(area)/sqrt(length(area)),3),3))%>%
  rename(Lifestage=lifestage, HPF=hpf)%>%
  kbl(caption="Descriptive statistics of larval size across ontogeny")%>%
  kable_classic(full_width=FALSE, html_font="Arial")%>%
  row_spec(0, bold = TRUE) 

Descriptive statistics of larval size across ontogeny

Lifestage	HPF	n	Mean	SE
Egg Fertilized	1	40	0.19	0.004
Embryo 1	5	40	0.184	0.003
Larvae 1	38	80	0.156	0.007
Larvae 2	65	40	0.158	0.003
Larvae 3	93	40	0.184	0.003
Larvae 4	163	40	0.193	0.004
Larvae 5	183	40	0.196	0.005
Larvae 6	231	40	0.312	0.014
Recruit 1	183	40	0.193	0.006
Recruit 2	231	33	0.27	0.021

#need to output to csv 
size%>%
  group_by(lifestage, hpf)%>%
  summarise(n=length(area),
            Mean=format(round(mean(area), 3), 3), 
            SE=format(round(sd(area)/sqrt(length(area)),3),3))%>%
  rename(Lifestage=lifestage, HPF=hpf)%>%
  write_csv(., "Mcap2020/Output/Physiology/larval_size_table.csv")

Plot data as a scatterplot

size$hpf<-as.factor(size$hpf)
size_plot<-size %>%
    ggplot(., aes(x = hpf, y = area)) +
    #geom_boxplot(outlier.size = 0) +
    geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
    geom_point(aes(fill=group, group=group), pch = 21, size=4, position = position_jitterdodge(0.1)) + 
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
    ylab(expression(bold(paste("Planar Size (mm"^2, ")"))))+
    ylim(0,1)+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=14),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); size_plot

#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30

Plot data as box plot

size_plot2<-size %>%
    ggplot(., aes(x = hpf, y = area)) +
    geom_boxplot(aes(color=group), outlier.size = 0, lwd=1) +
    geom_point(aes(fill=group), pch = 21, size=2, position = position_jitterdodge(0.1)) + 
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
    scale_color_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
    ylab(expression(bold(paste("Planar Size (mm"^2, ")"))))+
    ylim(0, 0.8)+
    theme_classic() + 
    #geom_text(label="A", x=1, y=2500, size=4, color="black")+ #egg
    #geom_text(label="A", x=2, y=2500, size=4, color="black")+ #embryo 1
    #geom_text(label="A", x=3, y=2500, size=4, color="black")+ #larvae 1
    #geom_text(label="AB", x=4, y=4100, size=4, color="black")+ #larvae 2
    #geom_text(label="AB", x=5, y=4100, size=4, color="black")+ #larvae 3
    #geom_text(label="AB", x=6, y=4100, size=4, color="black")+ #larvae 4
    #geom_text(label="BC", x=6.8, y=4500, size=4, color="black")+ #larvae 5
    #geom_text(label="CD", x=7.2, y=6500, size=4, color="black")+ #recruit 1
    #geom_text(label="D", x=7.8, y=6500, size=4, color="black")+ #larvae6
    #geom_text(label="D", x=8.2, y=8700, size=4, color="black")+ #recruit2
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=14),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14)
      ); size_plot2

Statistical analysis

Run lmer on cells per larvae by sampling point, specified by sequence of samples taken (life stage, hpf). Use tube ID as random effect.

size_model<-lmer(area~lifestage + (1|tube.ID), data=size)
summary(size_model)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: area ~ lifestage + (1 | tube.ID)
##    Data: size
## 
## REML criterion at convergence: -1223.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7656 -0.4461 -0.0426  0.3090  5.9137 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tube.ID  (Intercept) 0.000000 0.00000 
##  Residual             0.002973 0.05453 
## Number of obs: 433, groups:  tube.ID, 40
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          0.190150   0.008621 423.000000  22.056  < 2e-16 ***
## lifestageEmbryo 1   -0.005650   0.012192 423.000000  -0.463  0.64332    
## lifestageLarvae 1   -0.034300   0.010559 423.000000  -3.248  0.00125 ** 
## lifestageLarvae 2   -0.031825   0.012192 423.000000  -2.610  0.00937 ** 
## lifestageLarvae 3   -0.006350   0.012192 423.000000  -0.521  0.60277    
## lifestageLarvae 4    0.003025   0.012192 423.000000   0.248  0.80417    
## lifestageLarvae 5    0.006000   0.012192 423.000000   0.492  0.62290    
## lifestageLarvae 6    0.122175   0.012192 423.000000  10.021  < 2e-16 ***
## lifestageRecruit 1   0.003225   0.012192 423.000000   0.265  0.79152    
## lifestageRecruit 2   0.079395   0.012823 423.000000   6.192 1.41e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lfstE1 lfstL1 lfstL2 lfstL3 lfstL4 lfstL5 lfstL6 lfstR1
## lfstgEmbry1 -0.707                                                        
## lifestgLrv1 -0.816  0.577                                                 
## lifestgLrv2 -0.707  0.500  0.577                                          
## lifestgLrv3 -0.707  0.500  0.577  0.500                                   
## lifestgLrv4 -0.707  0.500  0.577  0.500  0.500                            
## lifestgLrv5 -0.707  0.500  0.577  0.500  0.500  0.500                     
## lifestgLrv6 -0.707  0.500  0.577  0.500  0.500  0.500  0.500              
## lifstgRcrt1 -0.707  0.500  0.577  0.500  0.500  0.500  0.500  0.500       
## lifstgRcrt2 -0.672  0.475  0.549  0.475  0.475  0.475  0.475  0.475  0.475
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

qqPlot(residuals(size_model))

## [1] 401 411

leveneTest(residuals(size_model)~lifestage, data=size)

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group   9  13.482 < 2.2e-16 ***
##       423                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(size_model, type=2)

## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## lifestage 0.91654 0.10184     9   423  34.253 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Violation in normality and variance assumptions. Conduct non-parametric test (Kruskal Wallis).

kruskal.test(area~lifestage, data=size)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  area by lifestage
## Kruskal-Wallis chi-squared = 156.39, df = 9, p-value < 2.2e-16

Significant difference in size between lifestages.

View posthoc comparisons for differences between lifestages.

emm = emmeans(size_model, ~ lifestage)
cld(emm, Letters=c(LETTERS)) #letter display

##  lifestage      emmean      SE   df lower.CL upper.CL .group
##  Larvae 1        0.156 0.00610  9.4    0.142    0.170  A    
##  Larvae 2        0.158 0.00862 38.0    0.141    0.176  AB   
##  Larvae 3        0.184 0.00862 38.0    0.166    0.201  AB   
##  Embryo 1        0.184 0.00862 38.0    0.167    0.202  AB   
##  Egg Fertilized  0.190 0.00862 38.0    0.173    0.208  AB   
##  Larvae 4        0.193 0.00862 38.0    0.176    0.211   B   
##  Recruit 1       0.193 0.00862 38.0    0.176    0.211   B   
##  Larvae 5        0.196 0.00862 38.0    0.179    0.214   B   
##  Recruit 2       0.270 0.00951 50.5    0.250    0.289    C  
##  Larvae 6        0.312 0.00862 38.0    0.295    0.330     D 
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 10 estimates 
## significance level used: alpha = 0.05 
## NOTE: Compact letter displays can be misleading
##       because they show NON-findings rather than findings.
##       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

pairs(emm)

##  contrast                   estimate     SE   df t.ratio p.value
##  Egg Fertilized - Embryo 1   0.00565 0.0122 38.0   0.463  1.0000
##  Egg Fertilized - Larvae 1   0.03430 0.0106 21.4   3.248  0.0859
##  Egg Fertilized - Larvae 2   0.03182 0.0122 38.0   2.610  0.2480
##  Egg Fertilized - Larvae 3   0.00635 0.0122 38.0   0.521  0.9999
##  Egg Fertilized - Larvae 4  -0.00302 0.0122 38.0  -0.248  1.0000
##  Egg Fertilized - Larvae 5  -0.00600 0.0122 38.0  -0.492  1.0000
##  Egg Fertilized - Larvae 6  -0.12218 0.0122 38.0 -10.021  <.0001
##  Egg Fertilized - Recruit 1 -0.00323 0.0122 38.0  -0.265  1.0000
##  Egg Fertilized - Recruit 2 -0.07940 0.0128 44.2  -6.186  <.0001
##  Embryo 1 - Larvae 1         0.02865 0.0106 21.4   2.713  0.2293
##  Embryo 1 - Larvae 2         0.02618 0.0122 38.0   2.147  0.5082
##  Embryo 1 - Larvae 3         0.00070 0.0122 38.0   0.057  1.0000
##  Embryo 1 - Larvae 4        -0.00868 0.0122 38.0  -0.712  0.9993
##  Embryo 1 - Larvae 5        -0.01165 0.0122 38.0  -0.956  0.9931
##  Embryo 1 - Larvae 6        -0.12782 0.0122 38.0 -10.484  <.0001
##  Embryo 1 - Recruit 1       -0.00887 0.0122 38.0  -0.728  0.9991
##  Embryo 1 - Recruit 2       -0.08505 0.0128 44.2  -6.626  <.0001
##  Larvae 1 - Larvae 2        -0.00248 0.0106 21.4  -0.234  1.0000
##  Larvae 1 - Larvae 3        -0.02795 0.0106 21.4  -2.647  0.2557
##  Larvae 1 - Larvae 4        -0.03732 0.0106 21.4  -3.535  0.0480
##  Larvae 1 - Larvae 5        -0.04030 0.0106 21.4  -3.817  0.0264
##  Larvae 1 - Larvae 6        -0.15648 0.0106 21.4 -14.819  <.0001
##  Larvae 1 - Recruit 1       -0.03753 0.0106 21.4  -3.554  0.0462
##  Larvae 1 - Recruit 2       -0.11370 0.0113 26.6 -10.067  <.0001
##  Larvae 2 - Larvae 3        -0.02548 0.0122 38.0  -2.089  0.5456
##  Larvae 2 - Larvae 4        -0.03485 0.0122 38.0  -2.858  0.1536
##  Larvae 2 - Larvae 5        -0.03782 0.0122 38.0  -3.102  0.0909
##  Larvae 2 - Larvae 6        -0.15400 0.0122 38.0 -12.631  <.0001
##  Larvae 2 - Recruit 1       -0.03505 0.0122 38.0  -2.875  0.1486
##  Larvae 2 - Recruit 2       -0.11122 0.0128 44.2  -8.666  <.0001
##  Larvae 3 - Larvae 4        -0.00937 0.0122 38.0  -0.769  0.9987
##  Larvae 3 - Larvae 5        -0.01235 0.0122 38.0  -1.013  0.9896
##  Larvae 3 - Larvae 6        -0.12853 0.0122 38.0 -10.541  <.0001
##  Larvae 3 - Recruit 1       -0.00958 0.0122 38.0  -0.785  0.9984
##  Larvae 3 - Recruit 2       -0.08575 0.0128 44.2  -6.681  <.0001
##  Larvae 4 - Larvae 5        -0.00298 0.0122 38.0  -0.244  1.0000
##  Larvae 4 - Larvae 6        -0.11915 0.0122 38.0  -9.772  <.0001
##  Larvae 4 - Recruit 1       -0.00020 0.0122 38.0  -0.016  1.0000
##  Larvae 4 - Recruit 2       -0.07637 0.0128 44.2  -5.950  <.0001
##  Larvae 5 - Larvae 6        -0.11618 0.0122 38.0  -9.528  <.0001
##  Larvae 5 - Recruit 1        0.00278 0.0122 38.0   0.228  1.0000
##  Larvae 5 - Recruit 2       -0.07340 0.0128 44.2  -5.719  <.0001
##  Larvae 6 - Recruit 1        0.11895 0.0122 38.0   9.756  <.0001
##  Larvae 6 - Recruit 2        0.04278 0.0128 44.2   3.333  0.0496
##  Recruit 1 - Recruit 2      -0.07617 0.0128 44.2  -5.935  <.0001
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 10 estimates

2. Symbiont Density

Symbiont density per individual

Import and manipulate data

# Cell count data
sym_counts <- read_csv("Mcap2020/Data/Physiology/CellDensity/symbiont.counts.csv")

sym_counts <- left_join(sym_counts, metadata)
sym_counts$hpf <- as.factor(sym_counts$hpf)

Calculate cells and normalize to either planar size (eggs through metamorphosed recruits) or surface area (attached recruits)

# Calculate mean counts for each sample
df <- sym_counts %>%
  dplyr::select(tube.ID, num.squares, matches("count[1-6]")) %>%
  gather("rep", "count", -tube.ID, -num.squares) %>%
  group_by(tube.ID, num.squares) %>%
  summarise(mean_count = mean(count, na.rm = TRUE))

#match in identifying information
df$lifestage<-sym_counts$lifestage[match(df$tube.ID, sym_counts$tube.ID)]
df$total.volume.ul<-sym_counts$total.volume.ul[match(df$tube.ID, sym_counts$tube.ID)]
df$num.individuals<-sym_counts$num.individuals[match(df$tube.ID, sym_counts$tube.ID)]
df$surface.area<-sym_counts$surface.area[match(df$tube.ID, sym_counts$tube.ID)]
df$hpf<-sym_counts$hpf[match(df$tube.ID, sym_counts$tube.ID)]
df$group<-sym_counts$group[match(df$tube.ID, sym_counts$tube.ID)]
df$lifestage<-as.factor(df$lifestage)
df$group<-as.factor(df$group)

# Normalize counts by homogenat volume and surface area
df <- df %>%
  mutate(cells.mL = mean_count * 10000 / num.squares,
         cells = cells.mL * (total.volume.ul/1000),
         cells.ind = cells / num.individuals, 
         cells.mm = cells / surface.area)

sym_counts<-df

Plotting

Plot data with mean and standard error for larvae through metamorphosis (these counts have cells/individual). Plot attached recruits separately, these values are in cells per mm2. We will plot cells per unit surface area for all stages in later analyses.

Display cells per individual.

sym_counts %>%
  filter(!group=="Attached Recruit")%>%
  ggplot(aes(x = lifestage, y = cells.ind, color = lifestage)) +
  labs(x = "",y = "Cell Density per larva") +
geom_jitter(width = 0.1) +                                            # Plot all points
  stat_summary(fun.data = mean_cl_normal, fun.args = list(mult = 1),    # Plot standard error
               geom = "errorbar", color = "black", width = 0.5) +
  stat_summary(fun = mean, geom = "point", color = "black") + # Plot mean
  theme_classic()

Display cell density per mm2 in attached recruit plugs. Plug 1 = 48 hps, Plug 2 = 72 hps, Plug 3 = 96 hps

sym_counts %>%
  filter(group=="Attached Recruit")%>%
  ggplot(aes(x = lifestage, y = cells.mm, color = lifestage)) +
  labs(x = "",y = "Cell Density per mm2") +
geom_jitter(width = 0.1) +                                            # Plot all points
  #stat_summary(fun.data = mean_cl_normal, fun.args = list(mult = 1),    # Plot standard error
               #geom = "errorbar", color = "black", width = 0.5) +
  #stat_summary(fun.y = mean, geom = "point", color = "black") + # Plot mean
  theme_classic()

Present means and standard error of each group and save summary table.

sym_counts%>%
  group_by(group, hpf, lifestage)%>%
  summarise(n=length(cells.ind),
            Mean=format(round(mean(cells.ind), 0), 0), 
            SE=format(round(sd(cells.ind)/sqrt(length(cells.ind)),0),0))%>%
  rename(Lifestage=group, HPF=hpf)%>%
  kbl(caption="Descriptive statistics of Symbiodiniaceae cell densities per larva across ontogeny")%>%
  kable_classic(full_width=FALSE, html_font="Arial")%>%
  row_spec(0, bold = TRUE) 

Descriptive statistics of Symbiodiniaceae cell densities per larva across ontogeny

Lifestage	HPF	lifestage	n	Mean	SE
Attached Recruit	183	Plug 1	3	NA	NA
Attached Recruit	231	Plug 2	2	NA	NA
Attached Recruit	255	Plug 3	3	NA	NA
Egg	1	Egg Fertilized	4	1472	125
Embryo	5	Embryo 1	4	1831	124
Embryo	38	Larvae 1	4	1371	242
Embryo	65	Larvae 2	4	2646	238
Larvae	93	Larvae 3	4	2692	347
Larvae	163	Larvae 4	4	2848	206
Larvae	183	Larvae 5	4	3474	134
Larvae	231	Larvae 6	4	5142	386
Metamorphosed Recruit	183	Recruit 1	4	4829	480
Metamorphosed Recruit	231	Recruit 2	4	6424	659

sym_counts%>%
  group_by(group, hpf, lifestage)%>%
  summarise(n=length(cells.mm),
            Mean=format(round(mean(cells.mm), 0), 0), 
            SE=format(round(sd(cells.mm)/sqrt(length(cells.mm)),0),0))%>%
  rename(Lifestage=group, HPF=hpf)%>%
  kbl(caption="Descriptive statistics of Symbiodiniaceae cell densities per mm2 across ontogeny")%>%
  kable_classic(full_width=FALSE, html_font="Arial")%>%
  row_spec(0, bold = TRUE) 

Descriptive statistics of Symbiodiniaceae cell densities per mm2 across ontogeny

Lifestage	HPF	lifestage	n	Mean	SE
Attached Recruit	183	Plug 1	3	21503	3084
Attached Recruit	231	Plug 2	2	30977	185
Attached Recruit	255	Plug 3	3	42014	3847
Egg	1	Egg Fertilized	4	NA	NA
Embryo	5	Embryo 1	4	NA	NA
Embryo	38	Larvae 1	4	NA	NA
Embryo	65	Larvae 2	4	NA	NA
Larvae	93	Larvae 3	4	NA	NA
Larvae	163	Larvae 4	4	NA	NA
Larvae	183	Larvae 5	4	NA	NA
Larvae	231	Larvae 6	4	NA	NA
Metamorphosed Recruit	183	Recruit 1	4	NA	NA
Metamorphosed Recruit	231	Recruit 2	4	NA	NA

#need to output to csv 
sym_counts%>%
  group_by(group, hpf, lifestage)%>%
  summarise(n=length(cells.ind),
            Mean=format(round(mean(cells.ind), 0), 0), 
            SE=format(round(sd(cells.ind)/sqrt(length(cells.ind)),0),0))%>%
  rename(group=group, HPF=hpf)%>%
  write_csv(., "Mcap2020/Output/Physiology/cell_density_table.csv")

Plot data as a scatterplot

sym_counts$hpf<-as.numeric(as.character(sym_counts$hpf))

symb_plot<-sym_counts %>%
    filter(!group=="Attached Recruit")%>%
    droplevels()%>%
    ggplot(., aes(x = hpf, y = cells.ind)) +
    #geom_boxplot(outlier.size = 0) +
    geom_smooth(method="lm", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
    geom_point(aes(fill=group, group=group), pch = 21, size=4, position = position_jitterdodge(5)) + 
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
    ylab(expression(bold(paste("Symbiont cells individual"^-1))))+
    ylim(0,9000)+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=14),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); symb_plot

#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30
#ATTACHED: #BA55D3

Plot data as box plot

symb_plot2<-sym_counts %>%
    filter(!group=="Attached Recruit")%>%
    droplevels()%>%
    ggplot(., aes(x = as.factor(hpf), y = cells.ind)) +
    geom_boxplot(aes(color=group), outlier.size = 0, lwd=1) +
    #geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
    geom_point(aes(fill=group), pch = 21, size=4, position = position_jitterdodge(0.2)) + 
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
    scale_color_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"), guide="none")+
    ylab(expression(bold(paste("Symbiont cells individual"^-1))))+
    ylim(0,9000)+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=14),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14)
      ); symb_plot2

Statistical analysis

Run ANOVA on cells per larvae by sampling point, specified by sequence of samples taken (life stage, hpf).

sym_ind_model_data<-sym_counts%>%
      filter(!group=="Attached Recruit")%>%
      droplevels()

sym_ind_model<-aov(cells.ind~lifestage, data=sym_ind_model_data)
summary(sym_ind_model)

##             Df    Sum Sq  Mean Sq F value   Pr(>F)    
## lifestage    9 102932981 11436998   25.06 1.34e-11 ***
## Residuals   30  13689458   456315                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

qqPlot(residuals(sym_ind_model))

## [1] 33 30

leveneTest(residuals(sym_ind_model)~lifestage, data=sym_ind_model_data)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  9  1.7683 0.1166
##       30

Both normality and homogeneity of variance pass.

There is a significant effect of lifestage on cell densities. View posthoc comparisons for differences between lifestages.

emm = emmeans(sym_ind_model, ~ lifestage)
cld(emm, Letters=c(LETTERS)) #letter display

##  lifestage      emmean  SE df lower.CL upper.CL .group
##  Larvae 1         1371 338 30      681     2061  A    
##  Egg Fertilized   1472 338 30      782     2161  A    
##  Embryo 1         1831 338 30     1141     2521  A    
##  Larvae 2         2646 338 30     1956     3335  AB   
##  Larvae 3         2692 338 30     2002     3382  AB   
##  Larvae 4         2848 338 30     2158     3538  AB   
##  Larvae 5         3474 338 30     2784     4163   BC  
##  Recruit 1        4829 338 30     4139     5519    CD 
##  Larvae 6         5142 338 30     4452     5831     D 
##  Recruit 2        6424 338 30     5734     7114     D 
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 10 estimates 
## significance level used: alpha = 0.05 
## NOTE: Compact letter displays can be misleading
##       because they show NON-findings rather than findings.
##       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

pairs(emm)

##  contrast                   estimate  SE df t.ratio p.value
##  Egg Fertilized - Embryo 1    -359.6 478 30  -0.753  0.9988
##  Egg Fertilized - Larvae 1     100.4 478 30   0.210  1.0000
##  Egg Fertilized - Larvae 2   -1174.0 478 30  -2.458  0.3298
##  Egg Fertilized - Larvae 3   -1220.4 478 30  -2.555  0.2816
##  Egg Fertilized - Larvae 4   -1376.3 478 30  -2.881  0.1554
##  Egg Fertilized - Larvae 5   -2002.0 478 30  -4.191  0.0073
##  Egg Fertilized - Larvae 6   -3669.9 478 30  -7.683  <.0001
##  Egg Fertilized - Recruit 1  -3357.6 478 30  -7.029  <.0001
##  Egg Fertilized - Recruit 2  -4952.2 478 30 -10.368  <.0001
##  Embryo 1 - Larvae 1           460.0 478 30   0.963  0.9923
##  Embryo 1 - Larvae 2          -814.4 478 30  -1.705  0.7842
##  Embryo 1 - Larvae 3          -860.8 478 30  -1.802  0.7288
##  Embryo 1 - Larvae 4         -1016.7 478 30  -2.128  0.5232
##  Embryo 1 - Larvae 5         -1642.4 478 30  -3.438  0.0469
##  Embryo 1 - Larvae 6         -3310.3 478 30  -6.930  <.0001
##  Embryo 1 - Recruit 1        -2998.0 478 30  -6.276  <.0001
##  Embryo 1 - Recruit 2        -4592.6 478 30  -9.615  <.0001
##  Larvae 1 - Larvae 2         -1274.4 478 30  -2.668  0.2316
##  Larvae 1 - Larvae 3         -1320.8 478 30  -2.765  0.1940
##  Larvae 1 - Larvae 4         -1476.7 478 30  -3.091  0.1013
##  Larvae 1 - Larvae 5         -2102.4 478 30  -4.401  0.0042
##  Larvae 1 - Larvae 6         -3770.3 478 30  -7.893  <.0001
##  Larvae 1 - Recruit 1        -3458.0 478 30  -7.239  <.0001
##  Larvae 1 - Recruit 2        -5052.6 478 30 -10.578  <.0001
##  Larvae 2 - Larvae 3           -46.4 478 30  -0.097  1.0000
##  Larvae 2 - Larvae 4          -202.3 478 30  -0.423  1.0000
##  Larvae 2 - Larvae 5          -828.0 478 30  -1.733  0.7685
##  Larvae 2 - Larvae 6         -2495.9 478 30  -5.225  0.0005
##  Larvae 2 - Recruit 1        -2183.6 478 30  -4.571  0.0027
##  Larvae 2 - Recruit 2        -3778.2 478 30  -7.910  <.0001
##  Larvae 3 - Larvae 4          -155.9 478 30  -0.326  1.0000
##  Larvae 3 - Larvae 5          -781.6 478 30  -1.636  0.8201
##  Larvae 3 - Larvae 6         -2449.5 478 30  -5.128  0.0006
##  Larvae 3 - Recruit 1        -2137.2 478 30  -4.474  0.0035
##  Larvae 3 - Recruit 2        -3731.8 478 30  -7.813  <.0001
##  Larvae 4 - Larvae 5          -625.8 478 30  -1.310  0.9435
##  Larvae 4 - Larvae 6         -2293.7 478 30  -4.802  0.0014
##  Larvae 4 - Recruit 1        -1981.3 478 30  -4.148  0.0082
##  Larvae 4 - Recruit 2        -3575.9 478 30  -7.486  <.0001
##  Larvae 5 - Larvae 6         -1667.9 478 30  -3.492  0.0415
##  Larvae 5 - Recruit 1        -1355.6 478 30  -2.838  0.1690
##  Larvae 5 - Recruit 2        -2950.2 478 30  -6.176  <.0001
##  Larvae 6 - Recruit 1          312.3 478 30   0.654  0.9996
##  Larvae 6 - Recruit 2        -1282.3 478 30  -2.684  0.2249
##  Recruit 1 - Recruit 2       -1594.6 478 30  -3.338  0.0590
## 
## P value adjustment: tukey method for comparing a family of 10 estimates

Output data to file.

sym_counts %>%
  write_csv(., file = "Mcap2020/Output/Physiology/calculated_densities.csv")

Symbiont density per unit size

Data manipulation and correlation

First, test for correlation between symbiont cell density and larval size to see if there is a relationship.

Generate data frame with summarised size and cell density information for each time point from eggs to metamorphosed recruits, because we have data for size and counts for each sample. We do not include attached recruits here yet, becuase we cannot calculate densities per individual.

#read in data frame generated in previous chunk 
sym_counts<-sym_counts%>%
  dplyr::select(tube.ID, lifestage, group, hpf, cells.ind, cells.mm)

#grab size data
area<-size%>%
  group_by(tube.ID)%>%
  summarise(mean_area=mean(area, na.rm=TRUE))

area$tube.ID<-as.factor(area$tube.ID)
sym_counts$hpf<-as.factor(sym_counts$hpf)

corr<-left_join(sym_counts, area)

Generate number of symbiont cells per mm^2 area for each tube.

corr<-corr%>%
  mutate(counts_area=cells.ind/mean_area)%>%
  mutate(counts_area=ifelse(is.na(counts_area), cells.mm, counts_area)) #add attached recruit data already calculated as cells per mm2

Plot correlation between cell counts (cells per individual) and size (area mm^2).

correlation<-corr %>%
    filter(!group=="Attached Recruit")%>%
    ggplot(., aes(x = mean_area, y = cells.ind)) +
    #geom_smooth(method="lm", se=TRUE, fullrange=TRUE, level=0.95, color="black", fill="gray") +
    geom_point(aes(fill=group), pch = 21, size=4) + 
    xlab(expression(bold(paste("Larval Size (mm"^2, ")")))) + 
    scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
    scale_color_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
    xlab(expression(bold(paste("Individual Size (mm"^2, ")"))))+
    ylab(expression(bold(paste("Symbiont cells individual"^-1))))+
    #ylim(0, 9000)+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=14),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14)
      ); correlation

Test relationship with a spearman correlation.

cor.test(corr$mean_area, corr$cells.ind, method=c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  corr$mean_area and corr$cells.ind
## S = 3704, p-value = 8.651e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6525328

Significant correlation between size and cell counts. r=0.37, p=0.017

Plotting

Plot cells per mm^2 as a boxplot.

#order for creating a legend for all plots 
corr$group <- factor(corr$group, levels = c("Egg", "Embryo", "Larvae", "Metamorphosed Recruit", "Attached Recruit"))

cells_size_plot<-corr %>%
    ggplot(., aes(x = hpf, y = counts_area)) +
    geom_boxplot(aes(color=group), outlier.size = 0, lwd=1) +
    geom_point(aes(fill=group), pch = 21, size=4, position = position_jitterdodge(0.4)) + 
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30", "#BA55D3"), guide="none")+
    scale_color_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30", "#BA55D3"))+
    ylab(expression(bold(paste("Symbiont cells mm"^-2))))+
    #ylim(2000, 35000)+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=14),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14)
      ); cells_size_plot

Plot as linear relationship.

cells_size_plot2<-corr %>%
    ggplot(., aes(x = as.numeric(as.character(hpf)), y = counts_area)) +
    #geom_boxplot(aes(color=group), outlier.size = 0, lwd=1) +
    geom_point(aes(fill=group, group=group), pch = 21, size=4, position = position_jitterdodge(0.4)) + 
    geom_smooth(method="lm", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30", "#BA55D3"), guide="none")+
    scale_color_manual(name="Lifestage", values=c( "#8C510A", "#DFC27D","#80CDC1", "#003C30", "#BA55D3"))+
    ylab(expression(bold(paste("Symbiont cells mm"^-2))))+
    #ylim(2000, 35000)+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=14),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14)
      ); cells_size_plot2

Analyze differences in normalized cell counts by timepoint.

model<-corr%>%
  #filter(!group=="Attached Recruit")%>%
  #droplevels()%>%
  aov(counts_area~lifestage, data=.)

qqPlot(residuals(model))

## [1] 41 30

corr%>%
  #filter(!group=="Attached Recruit")%>%
  #droplevels()%>%
  leveneTest(residuals(model)~lifestage, data=.)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group 12  1.1717 0.3399
##       35

summary(model)

##             Df    Sum Sq   Mean Sq F value   Pr(>F)    
## lifestage   12 3.543e+09 295228137   24.02 2.22e-13 ***
## Residuals   35 4.303e+08  12293680                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

View posthoc differences.

emm = emmeans(model, ~ lifestage)
cld(emm, Letters=c(LETTERS)) #letter display

##  lifestage      emmean   SE df lower.CL upper.CL .group  
##  Egg Fertilized   7746 1753 35     4187    11305  A      
##  Larvae 1         8720 1753 35     5161    12279  AB     
##  Embryo 1         9923 1753 35     6364    13483  ABC    
##  Larvae 4        14732 1753 35    11173    18291  ABCD   
##  Larvae 3        14756 1753 35    11197    18315  ABCD   
##  Larvae 6        16510 1753 35    12951    20069  ABCDE  
##  Larvae 2        16867 1753 35    13308    20426   BCDE  
##  Larvae 5        17732 1753 35    14173    21291    CDE  
##  Plug 1          21503 2024 35    17393    25613     DEF 
##  Recruit 2       23301 1753 35    19742    26860     DEF 
##  Recruit 1       25137 1753 35    21578    28696      EF 
##  Plug 2          30977 2479 35    25944    36010       FG
##  Plug 3          42014 2024 35    37905    46124        G
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 13 estimates 
## significance level used: alpha = 0.05 
## NOTE: Compact letter displays can be misleading
##       because they show NON-findings rather than findings.
##       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

pairs(emm)

##  contrast                   estimate   SE df t.ratio p.value
##  Egg Fertilized - Embryo 1   -2177.8 2479 35  -0.878  0.9994
##  Egg Fertilized - Larvae 1    -974.3 2479 35  -0.393  1.0000
##  Egg Fertilized - Larvae 2   -9120.9 2479 35  -3.679  0.0365
##  Egg Fertilized - Larvae 3   -7010.7 2479 35  -2.828  0.2348
##  Egg Fertilized - Larvae 4   -6986.0 2479 35  -2.818  0.2391
##  Egg Fertilized - Larvae 5   -9986.6 2479 35  -4.028  0.0149
##  Egg Fertilized - Larvae 6   -8764.6 2479 35  -3.535  0.0517
##  Egg Fertilized - Plug 1    -13757.4 2678 35  -5.137  0.0007
##  Egg Fertilized - Plug 2    -23231.6 3036 35  -7.651  <.0001
##  Egg Fertilized - Plug 3    -34268.5 2678 35 -12.797  <.0001
##  Egg Fertilized - Recruit 1 -17390.9 2479 35  -7.014  <.0001
##  Egg Fertilized - Recruit 2 -15555.2 2479 35  -6.274  <.0001
##  Embryo 1 - Larvae 1          1203.5 2479 35   0.485  1.0000
##  Embryo 1 - Larvae 2         -6943.1 2479 35  -2.800  0.2467
##  Embryo 1 - Larvae 3         -4832.9 2479 35  -1.949  0.7562
##  Embryo 1 - Larvae 4         -4808.2 2479 35  -1.939  0.7620
##  Embryo 1 - Larvae 5         -7808.8 2479 35  -3.150  0.1239
##  Embryo 1 - Larvae 6         -6586.8 2479 35  -2.657  0.3166
##  Embryo 1 - Plug 1          -11579.6 2678 35  -4.324  0.0067
##  Embryo 1 - Plug 2          -21053.8 3036 35  -6.934  <.0001
##  Embryo 1 - Plug 3          -32090.7 2678 35 -11.983  <.0001
##  Embryo 1 - Recruit 1       -15213.1 2479 35  -6.136  <.0001
##  Embryo 1 - Recruit 2       -13377.4 2479 35  -5.396  0.0003
##  Larvae 1 - Larvae 2         -8146.6 2479 35  -3.286  0.0921
##  Larvae 1 - Larvae 3         -6036.4 2479 35  -2.435  0.4450
##  Larvae 1 - Larvae 4         -6011.7 2479 35  -2.425  0.4512
##  Larvae 1 - Larvae 5         -9012.3 2479 35  -3.635  0.0406
##  Larvae 1 - Larvae 6         -7790.3 2479 35  -3.142  0.1259
##  Larvae 1 - Plug 1          -12783.0 2678 35  -4.773  0.0019
##  Larvae 1 - Plug 2          -22257.3 3036 35  -7.330  <.0001
##  Larvae 1 - Plug 3          -33294.2 2678 35 -12.433  <.0001
##  Larvae 1 - Recruit 1       -16416.6 2479 35  -6.622  <.0001
##  Larvae 1 - Recruit 2       -14580.9 2479 35  -5.881  0.0001
##  Larvae 2 - Larvae 3          2110.2 2479 35   0.851  0.9996
##  Larvae 2 - Larvae 4          2134.9 2479 35   0.861  0.9995
##  Larvae 2 - Larvae 5          -865.7 2479 35  -0.349  1.0000
##  Larvae 2 - Larvae 6           356.3 2479 35   0.144  1.0000
##  Larvae 2 - Plug 1           -4636.5 2678 35  -1.731  0.8683
##  Larvae 2 - Plug 2          -14110.7 3036 35  -4.647  0.0027
##  Larvae 2 - Plug 3          -25147.6 2678 35  -9.391  <.0001
##  Larvae 2 - Recruit 1        -8270.0 2479 35  -3.336  0.0824
##  Larvae 2 - Recruit 2        -6434.3 2479 35  -2.595  0.3499
##  Larvae 3 - Larvae 4            24.7 2479 35   0.010  1.0000
##  Larvae 3 - Larvae 5         -2975.9 2479 35  -1.200  0.9902
##  Larvae 3 - Larvae 6         -1753.9 2479 35  -0.707  0.9999
##  Larvae 3 - Plug 1           -6746.7 2678 35  -2.519  0.3935
##  Larvae 3 - Plug 2          -16220.9 3036 35  -5.342  0.0004
##  Larvae 3 - Plug 3          -27257.8 2678 35 -10.179  <.0001
##  Larvae 3 - Recruit 1       -10380.2 2479 35  -4.187  0.0098
##  Larvae 3 - Recruit 2        -8544.5 2479 35  -3.446  0.0638
##  Larvae 4 - Larvae 5         -3000.6 2479 35  -1.210  0.9896
##  Larvae 4 - Larvae 6         -1778.6 2479 35  -0.717  0.9999
##  Larvae 4 - Plug 1           -6771.3 2678 35  -2.529  0.3881
##  Larvae 4 - Plug 2          -16245.5 3036 35  -5.350  0.0004
##  Larvae 4 - Plug 3          -27282.5 2678 35 -10.188  <.0001
##  Larvae 4 - Recruit 1       -10404.9 2479 35  -4.197  0.0095
##  Larvae 4 - Recruit 2        -8569.2 2479 35  -3.456  0.0624
##  Larvae 5 - Larvae 6          1222.0 2479 35   0.493  1.0000
##  Larvae 5 - Plug 1           -3770.8 2678 35  -1.408  0.9658
##  Larvae 5 - Plug 2          -13245.0 3036 35  -4.362  0.0061
##  Larvae 5 - Plug 3          -24281.9 2678 35  -9.067  <.0001
##  Larvae 5 - Recruit 1        -7404.3 2479 35  -2.986  0.1733
##  Larvae 5 - Recruit 2        -5568.6 2479 35  -2.246  0.5671
##  Larvae 6 - Plug 1           -4992.8 2678 35  -1.864  0.8038
##  Larvae 6 - Plug 2          -14467.0 3036 35  -4.764  0.0020
##  Larvae 6 - Plug 3          -25503.9 2678 35  -9.524  <.0001
##  Larvae 6 - Recruit 1        -8626.3 2479 35  -3.479  0.0591
##  Larvae 6 - Recruit 2        -6790.6 2479 35  -2.739  0.2753
##  Plug 1 - Plug 2             -9474.2 3201 35  -2.960  0.1826
##  Plug 1 - Plug 3            -20511.1 2863 35  -7.165  <.0001
##  Plug 1 - Recruit 1          -3633.5 2678 35  -1.357  0.9741
##  Plug 1 - Recruit 2          -1797.8 2678 35  -0.671  1.0000
##  Plug 2 - Plug 3            -11036.9 3201 35  -3.448  0.0636
##  Plug 2 - Recruit 1           5840.7 3036 35   1.924  0.7712
##  Plug 2 - Recruit 2           7676.4 3036 35   2.528  0.3884
##  Plug 3 - Recruit 1          16877.6 2678 35   6.302  <.0001
##  Plug 3 - Recruit 2          18713.3 2678 35   6.988  <.0001
##  Recruit 1 - Recruit 2        1835.7 2479 35   0.740  0.9999
## 
## P value adjustment: tukey method for comparing a family of 13 estimates

Generate summary table of descriptive statistics.

#need to output to csv 
corr%>%
  group_by(group, hpf, lifestage)%>%
  summarise(n=length(counts_area),
            Mean_sym_mm2=format(round(mean(counts_area), 0), 0), 
            SE=format(round(sd(counts_area)/sqrt(length(counts_area)),0),0))%>%
  rename(Lifestage=group, HPF=hpf)%>%
  write_csv(., "Mcap2020/Output/Physiology/normalized_size_cells_summary.csv")

3. Correlations

Respirometry correlated with density

Load dataframes

resp<-read.csv("Mcap2020/Output/Respiration/mean_respiration.csv")
size<-read.csv("Mcap2020/Output/Physiology/larval_size_table.csv")
size_density<-read.csv("Mcap2020/Output/Physiology/normalized_size_cells_summary.csv")
density<-read.csv("Mcap2020/Output/Physiology/cell_density_table.csv")

Standardize column names.

resp<-resp%>%
    rename(Lifestage=group)%>%
    dplyr::select(Lifestage, Mean_R, Mean_P, Mean_GP, Mean_PR)

resp$Lifestage<-as.factor(resp$Lifestage)

levels(resp$Lifestage) <- list(`Larvae 2` = "Larvae2", `Larvae 3`  = "Larvae3", `Larvae 5`  = "Larvae5", `Larvae 6`  = "Larvae6", `Recruit 1`  = "Recruit1", `Recruit 2`  = "Recruit2")

size<-size%>%
  dplyr::select(Lifestage, Mean)%>%
  rename(Mean_Size=Mean)

size$Lifestage<-as.factor(size$Lifestage)

size_density<-size_density%>%
  dplyr::select(lifestage, Mean_sym_mm2)%>%
  rename(Mean_Size_Density=Mean_sym_mm2, Lifestage=lifestage)

size_density$Lifestage<-as.factor(size_density$Lifestage)

density<-density%>%
    rename(Lifestage=lifestage, Mean_Density=Mean)%>%
  dplyr::select(Lifestage, Mean_Density)

density$Lifestage<-as.factor(density$Lifestage)

Combine dataframes.

master<-full_join(resp, size, by=c("Lifestage"))

master<-full_join(master, size_density, by=c("Lifestage"))

master<-full_join(master, density, by=c("Lifestage"))

head(master)

##   Lifestage      Mean_R      Mean_P    Mean_GP  Mean_PR Mean_Size
## 1  Larvae 2 -0.01886538 0.011875684 0.03074106 1.767823     0.158
## 2  Larvae 3 -0.01364444 0.013637833 0.02728227 2.139497     0.184
## 3  Larvae 5 -0.01967560 0.032156372 0.05183197 2.705802     0.196
## 4  Larvae 6 -0.02208091 0.032023557 0.05410447 2.446370     0.312
## 5 Recruit 1 -0.05541750 0.002830864 0.05824836 1.120693     0.193
## 6 Recruit 2 -0.04029479 0.014010579 0.05430537 1.358576     0.270
##   Mean_Size_Density Mean_Density
## 1             16867         2646
## 2             14756         2692
## 3             17732         3474
## 4             16510         5142
## 5             25137         4829
## 6             23301         6424

Generate correlation matrix

Check for correlations between all variables.

pdf("Mcap2020/Figures/Physiology/Correlations.pdf", height = 9, width = 9)
g <-ggpairs(master, columns=2:8)
print(g)
dev.off()

## quartz_off_screen 
##                 2

There is a significant correlation between size normalized symbiont cell density and respiration. Plot this relationship.

Run Pearson correlation between respiration and size normalized cell densities.

cor<-master%>%
  drop_na()

cor.test(cor$Mean_Size_Density, cor$Mean_R, method=c("pearson"))

## 
##  Pearson's product-moment correlation
## 
## data:  cor$Mean_Size_Density and cor$Mean_R
## t = -8.9183, df = 4, p-value = 0.0008739
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9974514 -0.7890338
## sample estimates:
##        cor 
## -0.9757648

There is a significant correlation r=-0.98, p=0.0009.

r_corr_plot<-master %>%
    drop_na()%>%
    ggplot(., aes(x = Mean_Size_Density, y = Mean_R)) +
    #geom_smooth(method="glm", se=FALSE, fullrange=TRUE, level=0.95, color="gray", fill="gray") +
    geom_point(aes(fill=Lifestage), pch = 21, size=4) + 
    xlab(expression(bold(paste("Symbiont cells mm"^-2)))) + 
    scale_fill_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#80CDC1", "#80CDC1", "#003C30", "#003C30"))+
    scale_color_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#80CDC1", "#80CDC1", "#003C30", "#003C30"))+
    ylab(expression(bold(paste("Respiration (nmol O2) individual"^-1))))+
    #ylim(0, 9000)+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=14),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14)
      ); r_corr_plot

4. Generate Figures

Generate physiology panel with all variables of interest.

# extract the legend from one of the plots
legend <- get_legend(
  # create some space to the left of the legend
  cells_size_plot + theme(legend.box.margin = margin(1,1,1,1))
)

#remove legends from plots  
size_plot2<-size_plot2+theme(legend.position="none")
symb_plot2<-symb_plot2+theme(legend.position="none")
r_corr_plot<-r_corr_plot+theme(legend.position="none")
cells_size_plot_l<-cells_size_plot+theme(legend.position="none")

#assemble plots
all_plots<-plot_grid(size_plot2, symb_plot2, cells_size_plot_l, r_corr_plot, labels = c("A", "B", "C", "D"), label_size=18, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,0.8), align="h")

all_plots_legend<-plot_grid(all_plots, legend, rel_widths = c(4, 0.5), ncol=2, nrow=1)

ggsave(file="Mcap2020/Figures/Physiology/Physiology_figure.pdf", all_plots_legend, dpi=300, width=24, height=6, units="in")
ggsave(file="Mcap2020/Figures/Physiology/Physiology_figure.png", all_plots_legend, dpi=300, width=24, height=6, units="in")

Early life history physiology